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Abstract.  We  developed  an  automated  system  that  registers  chest  CT 
scans  temporally.  Our  registration  method  matches  corresponding  ana¬ 
tomical  landmarks  to  obtain  initial  registration  parameters.  The  initial 
point-to-point  registration  is  then  generalized  to  an  iterative  surface-to- 
surface  registration  method.  Our  “goodness-of-fit”  measure  is  evaluated 
at  each  step  in  the  iterative  scheme  until  the  registration  performance  is 
sufficient.  We  applied  our  method  to  register  the  3D  lung  surfaces  of  11 
pairs  of  chest  CT  scans  and  report  promising  registration  performance.^ 


1  Introduction 

Chest  computed  tomography  (CT)  has  become  a  well-established  means  of  diag¬ 
nosing  primary  lung  cancer  and  pulmonary  metastases  and  evaluating  response 
of  these  malignant  lesions  to  treatments.  Diagnosis  and  prognosis  of  cancer  gen¬ 
erally  depend  upon  growth  assessment  of  pulmonary  lesions  on  repeat  CT  stud¬ 
ies  [44,  29] .  Chest  CT  scans  can  also  be  used  for  lung  cancer  screening,  which  has 
been  proposed  by  some  but  is  still  controversial.  Bronchogenic  cancer  remains 
the  leading  cause  of  cancer  death  in  the  United  States,  killing  160,000  people 
a  year.  The  overall  5-year  survival  rate  is  now  15%  [20],  but  early  detection 
and  resection  can  improve  the  prognosis  significantly.  For  example,  the  5-year 
survival  rate  for  Stage  I  cancer  is  67%  [28]. 

Our  long-term  objective  is  to  develop  an  image  analysis  system  that  supports 
the  radiologist  in  detecting  and  comparing  pulmonary  nodules  in  repeated  CT 
studies  in  a  clinical  setting.  Such  a  system  must  solve  the  classical  problems 
in  medical  image  analysis  -  segmentation,  detection,  and  registration  -  for  the 
important  domain  of  chest  CT  images.  References  [5]  and  [19]  describe  our  pre¬ 
liminary  system  that  automatically  segments  the  thorax,  lungs,  and  structures 
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within  the  lungs,  and  detects  nodules  in  axial  chest  CT  images.  In  this  system, 
human  intervention  is  needed  to  register  the  CT  studies.  In  the  current  paper,  we 
focus  on  automating  the  registration  task  and  introduce  a  method  for  automatic 
three-dimensional  (3D)  alignment  of  lung  surfaces  in  repeat  CT  scans. 

A  large  body  of  literature  has  been  published  on  registration  techniques.  Ref¬ 
erences  [8]  and  [40]  provide  excellent  surveys.  Registration  methods  in  the  med¬ 
ical  image  domain  focus  primarily  on  the  brain,  e.g.,  [1,  9, 10, 12, 14,  21-23,  39, 
41,46],  but  also  other  organ  systems,  such  as  spine  [13,31],  foot  [34],  breast  [36, 
43],  and  prostate  [3]. 

Registration  of  images  within  the  same  modality,  e.g.,  CT  with  CT  or  mag¬ 
netic  resonance  (MR)  with  MR,  and  across  modalities  has  been  addressed.  CT 
images  have  been  correlated  with  MR  images  [22,  23,  33,  35]  and  positron  emis¬ 
sion  tomography  (PET)  with  MR  images  [18,23,30,32,33,37,38]  for  the  brain. 
MR  and  CT  scans  of  the  head  have  been  registered  to  the  patient  skin  surface 
depth  data  obtained  by  a  laser  range  scanner  [12].  Bone  scans  have  been  reg¬ 
istered  with  bone  films  [34]  for  the  foot.  CT  studies  of  the  thorax  have  been 
registered  to  PET  studies  [45]  and  bronchoscopic  images  [7]. 

Some  registration  approaches  ensure  that  data  sets  are  obtained  with  prospec¬ 
tive  attention  to  patient  positioning,  for  example,  through  head  fixation;  others 
retrospectively  reorient  image  sets  using  fixed  external  skin-surface  or  bone- 
implanted  fiducial  markers  [11,24,26].  Techniques  have  been  developed  to  reg¬ 
ister  points  to  points  [16]  or  surfaces  [4,25,26,30],  and  correlate  subimages  [2, 
17,42].  Often  only  a  small  misalignment  of  the  images  is  assumed  [2].  Other 
registration  methods  require  some  manual  input  to  compensate  for  rotational 
and  translational  differences  between  two  studies  [18,30]. 

Registration  of  chest  radiographs  has  been  addressed  by  Kano  et  al.  [17].  To 
the  best  of  our  knowledge,  an  automated  system  to  register  chest  CT  images 
temporally  has  not  been  developed  yet.  Registration  of  thoracic  CT  studies  is 
challenging,  since  patient  position  and  orientation  varies  each  time  a  study  is 
obtained.  Other  obstacles  are  differences  in  inspiratory  volumes.  The  patient’s 
thorax  is  imaged  when  the  patient  is  in  maximal  inspiration  for  the  entire  scan. 
Not  all  patients,  however,  can  comply  with  this  request. 

In  this  paper,  the  lung  surfaces  of  two  CT  scans  are  segmented  and  registered 
for  11  patients.  The  scans  were  taken  solely  for  clinical  reasons  and  without 
external  fiducial  markers  or  any  particular  attention  to  patient  position.  We 
first  describe  an  anatomical  landmark-based  registration  method  in  Sections  2.1 
-  2.3.  We  then  generalize  it  to  surface-to-surface  registration  in  Section  2.4.  In 
Section  2.5,  we  improve  the  method  using  an  iterative  registration  algorithm  that 
is  based  a  registration  scheme  by  Besl  [4].  We  then  report  registration  results 
for  11  pairs  of  chest  CT  scans  in  Section  3  and  conclude  with  a  discussion  in 
Section  4. 
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2  Methods 

Registration  techniques  relate  points  in  two  different  data  sets  to  each  other. 
The  geometric  nature  of  a  registration  transformation  can  be  described  by  its 
“elasticity”  [40] .  The  degree  of  elasticity  increases  from  rigid  to  affine  to  curved 
mappings.  Rigid-body  transformations  describe  translations  and  rotations  within 
the  image  plane.  Affine  transformations  map  straight  lines  to  straight  lines  and 
therefore  allow  scaling  and  skewing  of  the  data  in  addition  to  rotation  and  trans¬ 
lation.  Curved  transformations  map  straight  lines  into  curves,  for  example  with 
polynomial  mappings  or  other  deformations. 

In  this  section,  we  first  describe  a  registration  method  for  anatomical  land¬ 
marks  that  is  based  on  an  affine  point-to-point  transformation  that  models  3D 
rotation  and  translation,  as  well  as  scaling.  We  then  generalize  the  affine  transfor¬ 
mation  to  curve-to-curve  and  surface-to-surface  registration  and  finally  propose 
an  iterative  scheme  to  improve  the  surface-to-surface  registration. 

2.1  Registration  of  Anatomical  Landmarks 

Registration  techniques  determine  the  absolute  orientation  of  one  data  set  with 
respect  to  the  other.  When  used  in  the  computer  vision  or  photogrammetry  liter¬ 
ature,  the  term  “absolute  orientation”  generally  implies  that  the  3D  coordinates 
of  corresponding  points  in  the  two  different  data  sets  are  known  [15].  For  our 
3D  data  sets,  it  is  difficult  to  establish  the  anatomical  correspondence  of  most 
voxels,  even  for  a  human  observer.  We  therefore  use  the  voxels  that  make  up 
anatomical  landmarks,  such  as  spine  and  sternum,  for  our  initial  registration 
method.  Note  that  we  do  not  use  external  fiduciary  markers,  since  they  would 
be  impractical  in  a  clinical  setting. 

Bones  are  rigid  anatomical  features  that  can  be  registered  reliably  with  a 
3D  affine  transformation  that  models  rotation,  translation,  and  scaling.  In  par¬ 
ticular,  the  sternum  and  vertebrae  are  excellent  anatomical  landmarks,  because 
their  positions  are  relatively  fixed  within  the  chest.  We  estimate  the  3D  position 
of  anatomical  landmarks  by  identifying  the  landmarks  in  the  axial  images  and 
calculating  their  centroids. 

We  also  use  the  trachea  as  an  anatomical  landmark.  Although  trachea  po¬ 
sition  and  shape  change  with  respiration,  we  found  that  the  tracheal  centroids 
in  the  axial  images  serve  as  reliable  landmarks  for  registration  of  our  data  sets. 
Finally,  we  also  tested  the  use  of  structures  within  the  lungs,  for  example  nod¬ 
ules,  for  registration.  Figure  1  shows  how  the  centroids  of  sternum,  trachea,  and 
a  nodule  in  the  left  lung  are  registered  in  corresponding  axial  images  of  two  CT 
data  sets. 

We  describe  a  method  to  detect  anatomical  landmarks  in  Section  2.2  and  the 
3D  affine  landmark-to-landmark  transformation  in  Section  2.3. 

2.2  Correlation-Based  Recognition  of  Anatomical  Features. 

We  use  template  images  of  anatomical  landmarks,  for  example,  sternum,  verte¬ 
bra,  and  trachea,  as  shown  in  Fig.  1,  to  detect  these  landmarks  in  our  test  data. 
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Fig.  1.  On  the  left,  template  images  of  the  sternum,  trachea,  and  spine  are  shown.  In 
the  middle,  two  corresponding  CT  images  are  shown  that  are  part  of  two  different  CT 
studies  of  the  same  patient.  The  landmarks  chosen  for  the  initial  registration  are  the 
centroids  of  the  sternum,  trachea,  and  a  nodule  in  the  left  lung.  The  green  test  points 
in  study  2  must  be  matched  to  the  red  model  points  in  study  1.  On  the  right,  the  best 
affine  transformation  of  the  green  test  points  is  shown  in  blue. 


The  template  images  are  created  offline  by  manually  cropping  subimages  of  the 
features  out  of  a  training  data  set.  Although  the  features  look  slightly  different 
in  the  test  data,  training  and  test  data  generally  match  well.  This  is  particularly 
true  if  we  use  a  deformable  template  that  can  be  scaled  or  rotated. 

Let  a  describe  the  affine  parameters  position,  scale,  and  rotation  of  the  tem¬ 
plate.  We  use  the  normalized  correlation  coefficient  to  find  the  best  estimate  of 
the  affine  parameters.  In  our  previous  work  [6],  we  showed  that  the  statistically 
optimal  estimator  for  the  affine  parameters  takes  the  form  of  the  normalized  cor¬ 
relation  coefficient.  It  quantifies  how  well  the  measured  data  in  subimage  Iq{x,  y) 
matches  the  template  feature  in  q{x,  y;  a)  and  serves  as  a  match  measure  that  is 
“information  conserving”  because  it  exploits  all  the  measured  data  relevant  to 
the  feature’s  recognition.  The  normalized  correlation  coefficient  is  defined  by 

r(a)  =  — — — (A(a)  V  Iq{x,y)q{x,y;a)  -  mi{a)mg{a)),  (1) 

where  TO/(a)  =  '^Iq{x,y)  and  mq{a)  =  '^q{x,y;a)  are  the  respective  local  im¬ 
age  means,  erf  (a)  =  A{a)  Iq{x,  yf  -  (E  y)f  and  a‘^{a)  =  A(a)  E  q{x,  y, 
a)^  —  (E  9(2^1  y;  are  the  respective  local  variances,  and  where  the  sums  are 
computed  over  a  region  O  that  is  the  union  of  all  pixels  that  contain  the  expected 
feature  and  A  =  \0\  is  the  number  of  pixels  in  O. 

We  propose  a  hierarchical  method  of  search  for  the  global  peak  that  computes 
the  correlation  coefficient  at  different  resolution  scales.  For  position  estimation, 
for  example,  we  choose  a  grid-based  approach  which  samples  the  ambiguity  sur¬ 
face  at  every  10th  pixel,  then  around  local  peaks  at  every  5th  pixel,  and  then 
eventually  in  the  surrounding  of  maximum  peak  at  every  pixel.  The  ambigu¬ 
ity  surfaces  for  the  position  estimates  of  anatomical  features  have  global  peaks 
with  correlations  of  at  least  0.8,  which  lie  far  above  the  expected  correlation 
if  [r(a)]  =  0.  In  addition,  once  a  feature,  such  as  the  trachea,  is  found  in  an  axial 
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image,  the  search  space  for  the  same  feature  in  subsequent  images  can  be  reduced 
significantly.  In  addition,  the  template  feature  q  can  automatically  be  updated 
online  with  the  cropped  image  of  the  detected  feature  in  the  previous  slice.  This 
results  in  high  correlations  (above  0.9)  and  reliable  estimates  of  feature  position 
and  size. 

2.3  Three-Dimensional  AfRne  Point-to-Point  Registration. 

Given  a  voxel  x  in  study  1  and  a  voxel  p  in  study  2,  the  general  3D  affine 
transformation 

x  =  Ap-|-xo  (2) 

maps  p  into  x,  where  the  3x3  matrix  A  can  be  expressed  in  terms  of  nine 
parameters,  three  for  rotation,  three  for  scaling,  and  three  for  skewing.  Vector  xq 
describes  the  3D  translation  between  x  and  p. 

In  our  application,  the  orientation  of  the  patient’s  body  on  the  CT  table  is 
modeled  by  rotation  around  the  x,  y,  and  2  axes.  Parameters  that  model  scaling 
in  the  x-  and  y-dimensions  are  needed  if  the  field-of-view,  i.e.,  the  pixel-width- 
to-millimeter  ratio,  differs  between  two  studies.  Scaling  is  usually  uniform  in 
the  X-  and  y-  but  not  in  the  ^-dimension.  Scaling  in  z  is  due  to  the  differing 
slice  thickness  and  number  of  slices  in  the  two  studies.  Note  that  the  scaling 
parameters  are  determined  before  scan  acquisition.  We  therefore  do  not  need 
to  invert  for  the  scaling  parameters,  but  instead  directly  use  the  field-of-view 
and  collimation  information  provided  with  the  scan  data.  We  assume  that  the 
CT  scanner  does  not  introduce  skewing  and  preserves  the  Cartesian  (or  rectan¬ 
gular)  coordinates  of  3D  points.  Then  the  problem  of  finding  the  general  affine 
transformation  between  two  CT  studies  reduces  to  the  problem  of  finding  the 
rigid-body  transformation  between  the  two  studies  after  they  have  been  adjusted 
for  scaling  differences. 

The  rigid-body  transformation  T  maps  p  into  x, 

X  =  T(p)  =  Rp-fxo,  (3) 

where  the  orthonormal  3x3  matrix  R  rotates  p  into  vector  Rp,  which  is  then 
shifted  into  x  by  translation  vector  xq.  We  have  12  unknowns  (9  matrix  coef¬ 
ficients  and  3  translation  parameters)  and  only  3  linear  equations.  So  we  need 
at  least  4  corresponding  points  to  compute  the  unknown  transformation  pa¬ 
rameters.  If  we  impose  the  orthonormality  condition,  we  obtain  an  additional 
equation  and  therefore  only  need  3  corresponding  points.  Note  that  these  three 
points  must  not  be  collinear. 

Since  there  may  be  errors  in  the  measurement  of  the  points  or  in  the  corre¬ 
sponding  landmark  detection  algorithm,  a  greater  accuracy  in  determining  the 
transformation  parameters  can  be  obtained  if  more  than  three  points  are  used. 
Given  a  set  X  of  n  points  xi, . . .  ,x„  in  study  1  and  a  set  P  of  corresponding 
points  pi, . . .  ,p„  in  study  2,  our  goal  is  to  minimize  the  sum  of  square  residual 

^  l|eif  =  ^  ||x*  -  r(p,)f  =  ^  ||xi  -  Rp,  -  xof 

i—1 


errors 


(4) 
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with  respect  to  the  unknowns  R  and  xq.  A  closed- form  optimal  solution  to  this 
least-squares  problem  was  given  by  Horn  [16].  The  advantages  of  his  solution 
are  that  an  iterative  scheme  or  initial  guess  are  not  required.  The  best  possible 
transformation  is  computed  in  one  step.  An  additional  advantage  is  that  the 
solution  is  symmetric,  which  means  that  the  solution  that  transforms  P  into  X 
is  the  exact  inverse  of  the  solution  that  transforms  X  into  P. 

The  best  translation  vector  xq  is  the  difference  between  the  centroid  x  = 
l/nX]”=iXi  of  point  set  X  and  the  centroid  p  =  l/n-X^iLiPi  of  point  set  P 
rotated  by  rotation  R  : 

xo  =  x-R(p).  (5) 

Therefore,  the  translation  can  be  computed  easily  once  the  rotation  is  found. 
To  find  the  rotation,  the  coordinates  of  voxels  in  X  and  P  are  converted  into 
coordinates  of  voxels  in  X'  and  P'  of  coordinate  systems  that  are  originated 
at  the  respective  centroids,  e.g.,  x[  =  Xi  —  x  for  all  Xj  G  X.  This  reduces  the 
least-squares  problem  of  Eq.  4  to  a  minimization  of 


i—1  i—1  i—1  i—1 

with  respect  to  rotation  R  only,  or 


n 

max^xfRp'. 


The  solution  of  this  maximization  problem  is 

f  Qo  +  ^l-Qy-Qz  ‘^{QxQy  -  qoQz)  -  qoqy)  \ 

R  =  2{qyq^  -  qoqz,)  Qo  -  qI  +  -  Qz  "^iqyQz  -  qoqx) 

\  2{q^q^  -  qoqy)  2{q^qy  +  qoq^)  ql  -  ql  -  ql  +  ql ) 


(6) 


(7) 


(8) 


where  q=  {qo,qx,qy,qz)  the  unit  eigenvector  that  corresponds  to  the  maximum 
eigenvalue  of  the  symmetric  matrix 


/ 

^xx  T  ^yy  ^ zz 

^yz  ^zy 

^zx  ^xz 

^xy  ^yz 

\ 

^yz  ^zy 

^xx  ^yy  ^zz 

^xy  T  ^yz 

^zx  T  ^xz 

^zx  ^xz 

^xy  T  ^yz 

^xx  T  ^yy  ^zz 

^yz  T  ^zy 

V 

^xy  ^yx 

^zx  T  ^xz 

^yz  T  ^zy 

^xx  ^yy  ^zz 

/ 

and  Ski  is  the  kl-th  component  of  outer-product  matrix  S  = 

In  Figure  1,  the  centroids  of  the  sternum,  trachea,  and  a  nodule  in  study  1, 
shown  in  red,  make  up  the  model  set  X  =  {xi,X2,X3}  and  the  corresponding 
centroids  in  study  2,  shown  in  green,  make  up  the  test  set  P  =  {pi,  P2,  Ps}-  First 
we  compute  the  rotation  matrix  using  Eq.  8,  then  the  translation  parameters 
using  Eq.  5,  and  finally  the  transformation  T  of  P  into  X  using  Eq.  3,  which  is 
illustrated  in  blue  in  Fig.  1.  We  use  Eq.  4  to  report  the  “goodness  of  fit.” 
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2.4  Three-Dimensional  AfRne  Shape  Registration. 

In  this  paper,  we  focus  on  the  surface  of  the  lung  parenchyma  that  is  segmented 
from  the  full  3D  data  set  using  the  method  described  in  our  earlier  work  [19, 
5].  Figure  2  demonstrates  a  3D  view  of  the  lung  surface.  The  lung  parenchyma 
was  segmented  on  each  of  the  axial  images  of  a  low-dose  chest  CT  scan.  The 
scan  was  reconstructed  at  1.0  mm  intervals  and  with  a  1.25  mm  slice  thickness. 
Our  algorithm  identifies  the  lung  borders  of  the  parenchyma  and  uses  them  to 
construct  the  3D  surface  view  of  the  lung  parenchyma.  Our  goal  is  to  register  the 
lung  surfaces  segmented  on  an  initial  CT  scan  to  the  lung  surfaces  segmented 
on  another  CT  scan  of  the  same  patient  obtained  at  a  later  time. 


Fig.  2.  3D  visualization  of  the  lung.  Top  left  image:  coronal  view  of  both  lungs.  Top 
right  image:  top-down  axial  view  of  one  lung.  Bottom  images:  wire-frame  visualization 
of  one  lung. 


The  point-to-point  registration  algorithm  described  above  provides  the  ab¬ 
solute  orientation  of  one  point  set  P  with  respect  to  the  other  point  set  X.  It 
assumes  that  the  correspondence  between  points  and  Pi,  for  each  i,  has  been 
established.  For  certain  points,  for  example  the  centroids  of  the  sternum  in  cor¬ 
responding  axial  slices,  correspondence  can  be  determined  with  relatively  high 
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confidence,  but  for  other  point  pairs,  their  correspondences  are  not  as  easily  es¬ 
tablished.  For  example,  a  lung  border  point  in  the  apex  of  study  1  corresponds 
to  some  border  point  in  the  apex  of  study  2,  but  which  point  generally  cannot 
be  determined,  even  by  a  human  observer.  This  is  true  for  points  that  are  pixels 
on  curves  in  corresponding  2D  axial  images,  as  well  as  voxels  on  3D  surfaces. 
We  therefore  define  the  correspondence  C  of  two  points  on  different  curves  or 
surfaces  by  their  distances  to  each  other  and  the  rest  of  the  points.  In  particular, 
test  point  corresponds  to  model  point  =  C{pi),  if  their  Euclidean  distance 
is  the  shortest  among  all  distances  between  and  any  point  in  X,  i.e., 

C{pi)  =  Xj  for  which  ||xj  -  pj||  =  min  ||xfc  -  pi]].  (10) 

Note  that  C  is  not  a  symmetric  mapping,  i.e.,  the  corresponding  test  point  Pr  = 
C{xj)  of  model  point  xj  is  not  necessarily  p^,  since  the  shortest  distance  among 
all  distances  between  Xj  and  any  point  in  P, 

min  l|xj -p^ll,  (11) 

PaeP 

may  be  shorter  than  ||xj  —  pi|| .  Using  the  definition  for  correspondence  in  Eq.  10, 
we  can  match  two  curves  or  two  surfaces  to  each  other  that  contain  a  different 
number  of  pixels  or  voxels. 

The  correspondence  definition  in  Eq.  10  is  reliable  if  the  two  data  sets  are 
close  to  each  other,  in  particular,  if  they  have  been  registered.  This  creates  a 
paradoxical  situation:  we  would  like  to  register  corresponding  points,  but  need 
to  register  them  first  in  order  to  establish  their  correspondences.  To  resolve  this 
situation,  we  solve  the  registration  and  correspondence  problems  concurrently. 
We  developed  an  iterative  approach  based  on  Besl’s  iterative  closest-point  algo¬ 
rithm  [4]. 

We  first  detect  anatomical  landmarks  in  studies  1  and  2  and  compute  the  3D 
affine  transformation  that  registers  them  optimally,  as  described  in  Section  2.3. 
We  then  segment  the  lungs  [19,  5]  and  register  them  with  the  transformation  pa¬ 
rameters  computed  for  the  landmark  registration.  We  establish  correspondences 
by  computing  the  Euclidean  distances  between  all  point  pairs  of  the  two  data 
sets.  If  the  registration  error  is  too  large,  we  register  the  transformed  lung  bor¬ 
ders  in  study  2  to  the  lung  borders  in  study  1,  compute  the  new  correspondences 
and  registration  error,  and  then  iterate.  Once  the  error  is  sufficiently  small,  we 
terminate  the  process.  Convergence  of  this  iterative  algorithm  can  be  shown  [4] . 
The  pseudo-code  of  our  method  is  given  below. 

2.5  Registration  Code 

Function  LungRegistration  takes  as  inputs  3D  voxel  data  sets  CTstudyl  and 
CTstudyl  that  have  been  adjusted  for  field-of-view  differences,  and  a  parame¬ 
ter  threshold  that  is  used  to  decide  when  the  function  can  terminate  with  a 
sufficient  registration  result.  Function  LungRegistration  outputs  the  transfor¬ 
mation  parameters  for  translation  and  rotation.  Its  local  variables  are  defined 
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in  lines  2-6.  For  the  function  calls  in  lines  7  -  17,  we  use  C-style  notation 
to  distinguish  input  parameters,  for  example  lungl  in  line  14,  from  parameters 
that  change  during  the  function  call,  for  example  felungR,  fetranslation,  and 
ferotation  in  line  14. 

1  Function  LungRegistration  (CTstudyl,  CTstudy2,  threshold)  { 

2  float  error; 

3  voxelset  lungl,  lung2,  lungR; 

4  voxelset  landmarks 1,  landmarks2; 

5  voxel  translation; 

6  3x3matrix  rotation; 

7  DetectLandmarksC&landmarksl ,  &lEtndmarks2) ; 

8  RegisterLandmarks (landmarksl ,  lEtndmarks2, 

fetranslation,  ferotation) ; 

9  SegmentLungs (CTstudyl ,  felungl) ; 

10  SegmentLungs (CTstudy2,  &lung2) ; 

11  RegisterLungslnitially (lungl ,  lung2, 

translation,  rotation,  felungR) ; 

12  ComputeCorrespondencesAndError (lungl ,  lungR,  feerror) ; 

13  while  (error  >  threshold)  { 

14  RegisterLungs (lungl ,  felungR,  fetranslation,  ferotation) ; 

15  ComputeCorrespondencesAndError (lungl ,  lungR,  feerror); 

16  } 

17  OutputResults (translation,  rotation); 

18  } 

Variables  lungl  and  lung2  can  either  contain  the  voxels  that  make  up  the 
right  and  left  lungs  in  the  respective  CT  studies,  which  means  that  both  lungs 
will  be  registered  together,  or  they  can  contain  the  voxels  of  only  the  right 
or  the  left  lungs.  In  the  second  case,  LungRegistration  must  be  called  twice, 
once  for  the  right,  and  once  for  the  left  lung  registration.  Figure  3  shows  the 
results  of  joint  and  separate  lung  registrations  for  the  2D  case,  Fig.  4  the  results 
of  registering  3D  data  sets  that  contain  ten  left  lung  contours,  and  Fig.  5  the 
results  of  registering  the  full  3D  lung  surfaces  of  a  right  lung.  Visual  inspection 
shows  that  the  red  measured  and  blue  computed  points  match  well  and  confirms 
the  quantitative  error  analysis. 

3  Results 

Eleven  patients  with  cancer  diagnoses  and  pulmonary  nodules  were  selected,  who 
had  thoracic  CT  scans  for  clinical  indications  between  April  1993  and  August 
2000.  The  selected  patients  each  had  two  CT  studies,  and  a  total  of  22  CT  studies 
was  evaluated.  Ten  chest  CT  scans  had  been  performed  helically  on  GE  HiSpeed 
Advantage  machines  (GE  Medical  Systems,  Milwaukee,  WI)  according  to  the 
standard  departmental  protocol.  The  scans  were  obtained  from  the  lung  apices 
through  the  adrenal  glands  using  a  1:1  pitch  either  with  5mm  collimation  for  the 


10 


Betke,  Hong,  and  Ko 


Fig.  3.  Visualization  of  2D  lung  border  registration.  On  the  left,  two  corresponding  CT 
slices  are  shown  with  their  automatically  segmented  lung  contours  in  red  for  study  1  and 
green  for  study  2.  Three,  non-consecutive  iterations  during  the  registration  process  are 
illustrated  in  the  middle  and  right  images.  The  green  lung  borders  are  first  transformed 
to  the  dark  blue  contours,  then  from  there  to  the  light  blue  contours,  and  hnally,  after  a 
few  more  iterations,  to  the  white  contours.  In  the  middle  image,  the  resulting  contours 
are  shown  when  both  left  and  right  lungs  are  registered  together.  In  the  right  image, 
the  resulting  contours  are  obtained  from  two  registration  processes,  one  for  the  left 
lung  and  one  for  the  right.  Note  that  in  both  cases,  red  and  white  contours  match  well. 


Fig.  4.  Visualization  of  the  3D  registration  of  ten  left  lung  border  curves  with  views 
from  the  lung’s  side  (shown  in  the  image  on  the  top  left),  bottom  (bottom  left  image), 
and  top  (right  image) .  The  green  test  points  in  study  2  are  registered  to  the  red  model 
points  in  study  1  using  the  transformation  that  maps  green  points  to  blue  points.  The 
registration  error  is  minimal. 
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Fig.  5.  On  the  left,  wire  models  of  the  right  lung  for  study  1  are  shown  in  red  and 
study  2  in  green.  On  the  right,  the  transformed  lung  surface  of  study  2  is  shown  in 
blue.  It  matches  well  with  the  lung  surface  of  study  1,  shown  again  in  red. 


entire  study  or  10mm  collimation  with  5mm  collimation  through  the  hila.  One 
study  was  taken  on  a  multi-helical  Somatom  Volume  Zoom  CT  (Siemens  Medical 
Systems,  Iselin,  NJ).  The  scan  was  performed  using  a  1.0  mm  collimator  and 
reconstructed  in  1.25  mm  increments  at  1.0  mm  intervals.  Table  1  summarizes 
the  results. 


3.1  Registration  Speed 

The  advantage  of  the  initial  landmark  registration  is  that  it  significantly  increases 
the  speed  of  the  iterative  registration  process.  Note  that  we  could  have  merely 
guessed  an  initial  transformation  by  leaving  out  lines  7  and  8  and  switching 
lines  11  and  12  in  our  pseudo-code.  However,  depending  on  the  size  of  the  data 
sets,  registration  of  the  full  3D  lung  surfaces  without  initial  landmark  registration 
takes  in  the  order  of  hours  on  a  PC  with  a  866  MHz  Pentium  HI  processor. 
Instead  of  iterating  until  the  registration  error  falls  below  a  certain  threshold, 
we  fixed  the  number  of  iterations  to  be  25.  Figure  6  shows  how  the  error  decreases 
as  a  function  of  iteration  index.  For  one  of  our  data  sets,  after  25  iterations,  the 
registration  error  reduced  to  25%  of  the  initial  error.  For  a  data  set  reconstructed 
at  5/10/5  mm  thickness  creating  about  2  x  35  lung  slices,  25  iterations  take 
about  2  hours  on  average.  For  the  data  set  with  1.25  mm  thickness  creating 
about  2  X  200  lung  slices,  processing  25  iterations  took  more  than  3  days. 

Using  initial  landmark  registration,  a  much  smaller  number  of  iterations  is 
needed  to  produce  sufficient  registration  results.  Including  landmark  detection 
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Table  1:  Registration  Results 


Patient 

Months 

between 

Studies 

Reconstruc¬ 
tion  intervals 

in  mm 

No.  Lung 
Slices  in 
Studies  1,2 

Ratio  of 
Surface  Pts. 
Studies  2/1 

Rotation 
in  Euler 

Angles  (degrees) 

Trans¬ 

lation 

in  cm 

1 

2 

5/10/5 

31  31 

89% 

(0.4,  0.0,  6.3) 

2.9 

2 

5 

46  53 

114% 

(1.9,  -0.1,  3.9) 

3.3 

3 

4 

1 

196  202 

94% 

(1.0,  -4.7,  -0.5) 

1.6 

4 

74 

10/5/10,  5 

27  42 

154% 

(0.5,-0.9,  -9.0) 

2.6 

5 

2 

5 

42  49 

115% 

(0.2,  0.4,  -10.3) 

4.0 

6 

4 

10/5/10 

29  30 

88% 

(-0.4,  1.5,  -2.1) 

1.6 

7 

1 

•^2 

10/5/10 

29  30 

91% 

(0.4,  0.9,  -6.3) 

6.0 

8 

1 

5 

30  29 

92% 

(-1.5,  -1.2,  -2.9) 

3.7 

9 

4i 

10/5/10 

58  54 

98% 

(-0.5,  1.5,  -6.8) 

1.9 

10 

2 

5 

43  45 

105% 

^2.4,  -4.8,  12.8) 

9.5 

11 

5 

48  54 

78% 

^1.0,  1.6,  -12.2) 

7.2 

and  registration,  therefore  makes  our  algorithm  more  practical  and  cuts  the 
processing  time  significantly. 


Sum  of  Squared  Errors 


Iteration 


Fig.  6.  The  sum  of  squared  errors  is  shown  as  a  function  of  iteration  index  for  regis¬ 
tration  of  the  lung  surfaces  of  patient  2.  The  error  is  reported  in  units  of  10®  mm^. 


4  Discussion  and  Conclusions 

In  our  preliminary  system  [19,5],  image-to-image  registration  required  that  the 
lung  apices  were  identified  manually  on  the  two  studies.  Human  intervention  was 
also  needed  to  correlate  studies  with  different  collimation.  To  overcome  the  need 
of  manual  intervention,  we  developed  an  automatic  3D  registration  method  that 
matches  the  lung  surfaces  in  repeated  CT  studies. 
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Our  lung  registration  results  for  the  11  pairs  of  CT  scans  are  very  promising. 
The  performance  or  “goodness-of-fit”  of  our  registration  method  is  evaluated 
quantitatively  by  the  sum-of-squared-differences  measure  and  qualitatively  by 
visual  inspection.  In  the  future,  we  will  also  test  whether  our  system  can  reliably 
register  corresponding  nodules  in  repeated  chest  CT  scans  and  thus  become  a 
clinically  useful  tool  for  nodule  growth  assessment. 

We  presented  a  global  registration  method,  which  means  that  any  change  in 
a  transformation  parameter  influences  the  transformation  of  the  3D  data  set  as 
a  whole  [40] .  In  a  local  transformation,  such  a  change  influences  only  a  subset  of 
the  data.  In  the  future,  we  will  design  deformable  models  [27]  for  lung  surfaces 
in  order  to  model  local  transformations  that  are  due  to  differences  in  patient 
respiration.  We  will  use  the  deformable  model  parameters  that  register  lung 
border  surfaces  to  address  the  difficult  task  of  registering  structures  within  the 
lung.  This  will  require  modeling  the  deformable  shapes  of  nodules  in  3D  and 
also  modeling  nodule  position  as  a  function  of  lung  surface  deformation,  since 
nodules  may  move  within  the  lung  due  to  the  patient’s  respiration. 

Landmark  detection  and  registration  significantly  improve  the  speed  of  the 
registration  process.  Since  there  is  a  tradeoff  between  speed  and  precision  of 
registration,  we  will  test  the  impact  of  resolution  reduction  on  registration  per¬ 
formance,  in  particular,  nodule  registration.  We  will  also  investigate  whether 
initial  registration  of  a  larger  set  of  landmarks  will  improve  registration  preci¬ 
sion  and  speed. 

In  summary,  we  have  developed  a  3D  method  for  registration  of  lung  surfaces 
in  repeated  chest  CT  scans  and  applied  our  method  to  register  the  lungs  scans 
of  11  patients. 
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